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Abstract 

A transfer-matrix simulation scheme for the three-dimensional (d = 3) bond percolation is pre- 
sented. Our scheme is based on Novotny's transfer-matrix formalism, which enables us to consider 
arbitrary (integral) number of sites N constituting a unit of the transfer-matrix slice even for d = 3. 
Such an arbitrariness allows us to perform systematic finite-size-scaling analysis of the criticality 
at the percolation threshold. Diagonalizing the transfer matrix for iV~ = 4, 5, . . . , 10, we obtain an 
estimate for the correlation-length critical exponent u = 0.81(5). 
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I. INTRODUCTION 



The transfer-matrix method has an advantage over the Monte Carlo method in that 
it provides information free from the statistical (sampling) error and the problem of slow 
relaxation to thermal equilibrium. On one hand, the tractable system size is severely limited, 
because the transfer-matrix size increases exponentially as the system size enlarges. Such a 
difficulty could be even more serious for "geometrical" problems such as the percolation, for 
which the configuration space is much larger than that of the Ising model, for example. 

A transfer-matrix approach to the percolation in two dimensions (d = 2) was initiated by 
Derrida and Vannimenus Q]. They treated the system sizes (transfer-matrix strip widths) 
up to N = 5. Performing an extensive phenomenological renormalization group (finite-size 
scaling) analysis [2], they estimated the correlation-length critical exponent as v = 1.2 ~ 1.4; 
the variance is due to the choice of the boundary conditions. Their result is quite consistent 
with the exact value v = 4/3, indicating that the transfer-matrix approach to percolation 
would be promising. Because the transfer-matrix data are free from the statistical error, the 
data allow us to take its numerical derivative, which provides valuable information as to the 
subsequent finite-size-scaling analysis. 

It turned out, however, that its naive extension to the d = 3 case is rather problematic; 
we refer to Sec. 4.4 of Ref. ^ for an overview. Actually, for d — 3, as the system size (linear 
dimension) L enlarges, the number of constituent sites N(= L 2 ) of the transfer-matrix unit 
soon exceeds the limit of the available computer resources. 

The aim of this paper is to develop an improved version of the transfer-matrix formalism 
for the d = 3 bond percolation. For that purpose, we utilize Novotny's idea, which has been 
applied successfully to various Ising models in d < 7 □□□□□□ . His formalism stems on 
a very formal expression for the transfer-matrix elements. It enables us to consider arbitrary 
(integral) number of constituent sites even for d > 3. Owing to this arbitrariness, we are 
able to treat a variety of system sizes, and perform systematic finite-size scaling analysis of 
the numerical data. In this paper, we demonstrate that his idea is applicable to the d = 3 
bond percolation. 

The rest of this paper is organized as follows. In Sec. |Hj we formulate the transfer-matrix 
scheme for the d = 3 bond percolation. We place an emphasis how we adapted Novotny's 
idea to the bond-percolation problem. The simulation results are shown in Sec. IHII Taking 
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an advantage that we can treat various system sizes, we manage the phenomenological 
renormalization group 2] (finite-size scaling) analysis. Thereby, we obtain an estimate for 
the correlation-length critical exponent v = 0.81(5). In Sec. I1VI we present the summary 
and discussions. 



II. CONSTRUCTION OF THE TRANSFER MATRIX FOR THE THREE- 
DIMENSIONAL BOND PERCOLATION 

In this section, we set up the transfer-matrix formalism for the d = 3 bond percolation. 
As mentioned in Introduction, we adopt Novotny's idea 0. So far, his idea has been applied 
to various types of the Ising models in d < 7. Here, we show that his idea is also applicable 
to the bond percolation. At the end of this section, we argue a conceptual difference from 
the original Novotny method. 



A. Configuration space 

Above all, we need to set up the configuration space so as to represent the transfer- 
matrix elements explicitly. The bases of the configuration space should specify all possible 
connectivities among the N sites which constitute a unit of the transfer matrix, namely, 
a cross section of the transfer-matrix bar; see Fig. ^ In the figure, we also presented a 
drawing of an example of connectivity among the N = A sites. 

As shown in the figure, an integer index a = 1,2, ... ,N specifies the position of the 
constituent sites of a transfer-matrix unit. In order words, the transfer-matrix unit is of one- 
dimensional structure rather than two-dimensional one. Such a feature might be confusing, 
compared with the drawing in Fig. ^ (a), where the transfer- matrix unit is drawn as a 
rectangular shape with the edges y/~N. Actually, the dimensionality is lifted to d — 2 
afterward by introducing the A/iVth-neighbor long-range interactions among the N sites. 
We will explain this scheme explicitly in the next subsection. Here, for the time being, we 
consider that these N sites are arranged into a one-dimensional structure. 

In order to specify the connectivity among the N sites, we accepted the following matrix- 
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based representation; 



1 for a pair of connected sites (a, 0) 
otherwise 



(1) 



(2) 



The index i runs over all possible connectivities among the N sites. For example, the 
connectivity of Fig. d (b) is represented by, 

/ 1 1 1 \ 
110 1 
10 
\1 1 1/ 

Let us mention a few remarks: The authors in Ref. [jj, accepted more elaborated repre- 
sentation scheme. Namely, they specified whether a cluster is connected with the "origin" 
or that a cluster is an isolated one. Here, the "origin" stands for an edge of the transfer- 
matrix bar, An advantage of their extended representation is that one only needs to evaluate 
the largest eigenvalue of the transfer matrix to obtain the correlation length. Here, how- 
ever, we did not accept their representation scheme. Correspondingly, we calculated the 
sub-dominant eigenvalue together with the dominant one in order to calculate the corre- 
lation length. This task is not so computationally demanding, and it renders significant 
simplification of the formalism mentioned below. 



B. Explicit formula of the transfer-matrix elements for the d = 3 bond percolation 

In this subsection, we present the explicit formula for the transfer-matrix elements. We 
consider an anisotropic bond percolation on the cubic lattice. Namely, we set the percolation 
probabilities p», p±i andp^2 independently for the respective bond directions along the d = 3 
Cartesian axes. Correspondingly, we factorize the transfer matrix into the following three 
components; 

T = T ± (v,p± 2 )T ± (l,p ±1 )T l \(p ll ), (3) 

with, 

v = VN. (4) 

(We followed the idea of Derrida and Vannimenus, who decomposed the transfer-matrix 
for the d — 2 percolation into two factors jl].) The components Tji(p||), T±(l,p ±1 ) and 
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T±(v,p±2) denote the transition probabilities due to the longitudinal bond, intra-cluster 
nearest-neighbor bond, and the intra-cluster vth-neighbor bond, respectively. (Here, the 
"cluster" stands for the N sites constituting a unit of the transfer-matrix slice, and the 
"longitudinal" direction is parallel to the transfer-matrix bar; see Fig. [I]) 

In other words, the product of two components Tj|(p||)Tj_(l,pj_i), namely, with Tj_(t>,£>_|_ 2 ) 
ignored, should yield the transfer matrix for the d = 2 bond percolation. The remaining 
factor T±(v,p±2) lifts the dimensionality to d — 3. 

We present the explicit formulas for each component. First, for simplicity, we consider 
the longitudinal part Tji . (This is essentially the same as the horizontal factor Mh appearing 
in the formalism [if for the d = 2 percolation.) Our formula for the elements of Tji is given 
by, 

l T \\(P\\)]ij = ^2p({ j «},P\\) (at, m({ J a }) ® %) ■ (5) 

Here, the summation runs over all possible random-bond configurations {J a } with 

either J a = (unoccupied bond) or 1 (occupied bond) for a — 1, 2, . . . , N. The probability 
p({ J a},P\\) is given by, 

P({Ja},Pl\)=P^(l-p\\) N - Nj -\ (6) 

with the number of occupied bonds iVj =1 . The "random bond" matrix m({ J a }) is a diagonal 
N x N matrix with the diagonal elements [m({Ja})}i3f3 = J/3- The operation ® denotes the 
matrix product, 

N 

[a <g> b) a/3 = ^ ' a "7 A & 7/3> ( 7 ) 

7=1 

with the logical product A in the Boolean algebra. (The summation ^ gives 1, unless all 
the summands are zero.) The product (a i; aj) accounts for the orthogonality of the matrices; 
namely, 

(a,i,aj) = 8ij, (8) 

with Kronecker's symbol 5ij. 

As would be apparent from the above formula (0), the transfer-matrix element [Tji]ij 
stands for the transition probability from the initial configuration aj to the final configuration 
di through the longitudinal random-bond percolation {J a }. 

Second, we turn to considering the transverse component T±(w,p±). This component 
accounts for the intra-cluster wth-neighbor random-bond percolation with the percolation 
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probability p±. This factor is the most significant part in our formalism. We propose the 
following formula for T±(w,p±); 

[T±(w,p±)] tj = ^2p({J a },P±)tij(w,{J a }). (9) 
The transition amplitude is given by, 

N 

tij(w, {J a }) = ^ fp( W ) m p({ J *}) © a i) ■ ( 10 ) 

13=1 

Here, the symbol © denotes an operation, 

TV 

[a b\ a p = } j 'day V b ia , (11) 

7=1 

with the logical summation V in the Boolean algebra. The "/3th-neighbor-random-bond 
matrix" mp({J a }) is given by the formula, 

mp{{J a })=m{{J a })®sp, (12) 

with the shift operator, 

J 1 for 1 -5 = (3 mod N 
[sphs = < ■ (13) 

I otherwise 

The operation sp shifts the diagonal random-bond operator m({ J a }) to an off-diagonal one, 
which now represents the /3th-neighbor random bonds. That is, the operation mp({ J a }) (Bdj 
introduces new intra-cluster /3th-neighbor-random-bond percolation in adding to the initial 
connectivity a,j. 

In order to implement the wth-neighbor intra-cluster percolation with a non-integral value 
of w, we need to average over all sectors (3 = 1,2, . . . , N with an appropriate weight fp(w); 
see Eq. (JTUJ. We propose that the weight should be given by the u>th-order power of the 
shift operator; 

u(w) = [( ai rv (M) 

As would be apparent from the definition, the operator (si) w generates the translational 
shift of the distance w, and the factor fp(w) picks up the amplitude of each sector (3. 
Hence, the resulting formula, Eq. (jlOjl . should be the transition amplitude from aj to a, by 
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the wth-neighbor random-bond percolation. Hence, the product T±(v,p±2)T±(l,p±i) with 
v = yN introduces a two-dimensional intra-cluster percolation among the N sites. As noted 
previously, a crucial point is that there is no restriction to the number of constituent sites 
N. 

r-i 

Lastly, let us argue a conceptual difference from the original Novotny method |4| for 
the Ising ferromagnet. In the original method, the translation operator (si) w acts on the 
configuration space, which has huge dimensionality. On the contrary, in our formalism, 
the operator (si) w is a mere N x N matrix. Hence, from the viewpoint of the computer 
programing, the present formalism for the percolation is even simpler than the original 
method. (Here, the generation of the list of connectivity {a{\ is the most time-consuming.) 

III. NUMERICAL RESULTS 

In the above section, we set up the transfer-matrix formalism for the d = 3 bond percola- 
tion. In this section, we present the numerical results by means of the exact diagonalization 
of the transfer matrix. We consider the anisotropic bond percolation. The anisotropy pa- 
rameter R governs the mutual ratios of the percolation probabilities; namely, 

R- x P\\ = R- 1 px.x=Px.2=P- (15) 

We consider the anisotropy ratio R as a freely tunable parameter to stabilize the finite-size 
corrections. Diagonalizing the system sizes N = 4, 5, ... , 10, we analyze the percolation 
transition in terms of the phenomenological renormalization group [2] method. Note that 
the linear dimension of the system L is given by the relation, 

L = VN, (16) 

as shown in Fig. ^ 

A. Percolation threshold p c 

In Fig. |21 we plotted the scaled correlation length £/L for the percolation probability p 
with the fixed anisotropy parameter R = 3.3. (Afterward, we explain how we adjusted the 
anisotropy parameter to R = 3.3.) The symbols +, x, *, □, ■, o, and • denote the system 
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sizes N = 4, 5, 6, 7, 8, 9, and 10, respectively. The correlation length £ is calculated by the 
formula £ = l/ln(Ai/A2) with the dominant Ai and the sub-dominant A2 eigenvalues of the 
transfer matrix. 

The intersection point of the curves in Fig. |2] indicates the location of the critical point 
(percolation threshold) p c . (The scaled correlation length £/L is invariant with respect to 
the system size N at p = p c .) Hence, we observe that a percolation transition occurs around 
p » 0.096. 

In order to determine the critical point p c more precisely, in Fig. El we plotted the 
approximate transition point p c {Li, L 2 ) for (2/(Li + L2)) 2 . Here, the approximate transition 
point p c (Li,L 2 ) denotes the location of the intersection point of the curves for a pair of 
(Li, L 2 ). The symbols +, x, and * show that the differences of the system sizes are Ni—N 2 = 
1, 2, and 3, respectively; note that the relation Li 2 = v/ N lj2 holds. As indicated in Fig. 
we survey several values of the anisotropy parameter such as R = 2.8, 3.3 and 3.8. Thereby, 
we notice that the finite-size corrections to p c are suppressed on setting R = 3.3. The 
least-square fit to the data for R = 3.3 yields the critical point p c = 0.0958(27) in the 
thermodynamic limit L — > 00. 

Let us argue the role of the anisotropy parameter R. First of all, it is worthwhile that 
the system size along the transfer-matrix (longitudinal) direction is infinite, whereas the 
transverse system sizes are both finite L < vTO; see Fig. In this sense, it is by no means 
necessary to consider the isotropic condition p±\ = pi_i = p\\ specifically. Hence, we consider 
that the ratio R = p\\/p± 2 is a tunable parameter. Practically, we found that the finite-size 
corrections improve for large R. Second, we need to remedy the dimensionality d = 3 by 
adjusting the ratio of the intra-cluster interactions R = p±i/p± 2 . That is, the "effective 
dimension" P, |(| can deviate slightly from d — 3, at least, for small system sizes. (This 
deviation deteriorates the finite-size-scaling analysis.) Note that basically, the backbone 
structure of Novotny's transfer matrix is of d = 2, and the dimensionality is lifted to d = 3 
by introducing the long-range interactions among the intra-cluster sites. In other words, it 
is not quite obvious that the dimensionality d = 3 is realized precisely, at least, for small 
N. Hence, in order to remedy this dimensionality deviation, we should tune ^, 0] the intra- 
cluster-interaction ratio R = p±i/p± 2 . (Note that for large R, the component T±(l,p±x) 
dominates T^(v,p± 2 ), and the dimensionality reduces to d — 2. For a certain moderate 
value of R, the dimensionality would approach d — 3.) More specifically, we adjusted R 
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so as to improve the finite-size-scaling behaviors of p c (Li,L 2 ) p, |6| as shown in Fig. El In 
this respect, there might exist alternative parameterization schemes other than the present 
one. Here, however, we accepted the simplest parameterization scheme py = p±i = Rp±2 
involving a single tunable parameter R. 



B. Correlation- length critical exponent v 

In this subsection, we study the criticality at p — p c . In Fig. |U we plotted the approximate 
correlation-length critical exponent 



dp dp 



(17) 



p=p c 

for (2/(Li + L 2 )) 2 . Here, we set p c = 0.0958 and R = 3.3. The symbols +, x, and * show 
that the differences of the system sizes are Ni — N% — 1, 2, and 3, respectively. Note that this 
formula contains p derivative, which is readily calculated with the transfer-matrix method 
very precisely. (The transfer- matrix data are free from the statistical error.) 

We see that the data align rather satisfactorily. The least-square fit to these data yields 
an extrapolated value v = 0.812(15) to the thermodynamic limit L — > oo. Similarly, we 
obtain v = 0.813(21) by omitting the L — 10 data. Thereby, we confirm that the data are 
almost converged. 

In order to check the reliability of v further, we try to manage alternative extrapolation 
schemes: First, we replace the abscissa in Fig. |3]with the refined one (2/(Li + L2)) U)+1 ^ U 
jlp| |. where we used uj = 1.61(5) and v = 0.89(2) reported in Ref. Thereby, we arrive at 
v = 0.811(15). This result indicates the stability of v with respect to p c . Second, replacing 
the abscissa in Fig. H]with (2/(Li + £2))^, we obtain v = 0.853(19). Actually, this refined 
extrapolation yields an "improved" value for v. However, for the sake of self-consistency, we 
do not accept this refined extrapolation method, and consider it as a reference. Lastly, setting 
the values of the anisotropy parameter as R = 2.8 and R = 3.8, we obtain v = 0.853(27) 
and v = 0.771(11), respectively. These results again confirm the stability of v satisfactorily. 

Recollecting these results, we estimate the correlation-length critical exponent as, 

^ = 0.81(5), (18) 

with an expanded error margin. 
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Let us recollect a number of recent estimates for v determined with other approaches: 
From the series expansion method, Dunn et al. obtained v = 0.83(5). On the other hand, 
;he Monte-Carlo studies have reported v = 0.89(2) Q, v = 0.8765(16) 13|, v = 0.893(40) 
l^ . and v = 0.868(11) 15]. These estimates and ours are consistent with each other within 
the error margins. Nevertheless, we stress that our motivation is not necessarily directed to 
the accurate estimation of the critical indices. In the last section, we address an extended 
remark on the potential applicability of our scheme and future perspective. 

Lastly, in Fig. we present the approximate (3 function [l^ . 

0{Ll , l 2 ) = '-MfM/iM)/KW (19) 

\l °% l)a % i, li(L 1 )lw 

for R = 3.3. (Note that this formula also contains the derivatives, and it is hardly accessible 
by other approaches.) The symbols +, x, and * show that the pairs of system sizes are 
(Ni,N 2 ) = (6,8), (7,9), and (8,10), respectively. The beta function provides rich infor- 
mation on the overall feature of the criticality. The zero point fi(p)\ p=Pc = indicates the 
location of the transition point p c , and the slope at the transition point yields the inverse of 
v. In the figure, we presented a slope — {p — p c )/v with p c = 0.0958 and v = 0.81 determined 
above. The slope well describes the behavior of the beta function in the vicinity of p = p c . 
However, in a closer look, the beta function bends convexly, indicating that non-negligible 
corrections to scaling do exist. Possibly, such severe corrections are reflected in Fig. where 
we observe pronounced finite-size corrections to v. However, these corrections are fairly sys- 
tematic so that we could manage the extrapolation to L — > oo rather unambiguously. This 
is an advantage of Novotny's method, with which a variety of system sizes are available even 
for the case of d = 3. 



C. Isotropic case p» = p±i = p±2 

In Fig. El we plotted the beta function (3(L\, L 2 ) for R = 1 (isotropic case). The symbols 
+ , x, and * show that the pairs of system sizes are (Ni,N 2 ) = (6,8), (7,9), and (8, 10) 
respectively. We also presented a slope — (p — p c )/v with p c = 0.2488126 and v = 0.89 [ll 
as a dashed line. We see that our data and the slope behave similarly in the vicinity of 
p = p c . We obtain the correlation-length critical exponent v « 0.64 from a pair of iV = 8 
and 10. However, for small N, the data scatter, and eventually, even the zero point of 
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the beta function disappears. Because of this irregularity, we cannot manage systematic 
extrapolation to the thermodynamic limit. In this sense, the anisotropy parameter R is 
significant in order to stabilize the finite-size corrections of the transfer-matrix data as 
demonstrated in the preceding subsections. 

IV. SUMMARY AND DISCUSSIONS 

We developed a transfer-matrix formalism for the d = 3 bond percolation. Our formalism 
is based on Novotny's idea |4], which has been applied to the Ising models in high dimensions 
d < 7. We demonstrated that his idea is also applicable to the d = 3 bond percolation. A 
key ingredient of his method is that we can treat arbitrary number of sites constituting a 
unit of the transfer-matrix slice; see Fig. 

Diagonalizing the transfer matrix for the system sizes N = 4,5,. ..,10, we studied the 
criticality of the percolation transition. We found that the numerical data are well described 
by the finite-size-scaling theory, and thereby, we obtained an estimate for the correlation- 
length critical exponent v = 0.81(5). Here, we tuned the anisotropy parameter to R = 3.3 
in order to reduce the finite-size corrections. Because the system size along the transfer- 
matrix direction is infinite, it is by no means necessary to consider the isotropic limit R = 1 
specifically. 

The aim of this paper is to demonstrate an applicability of the transfer-matrix method 
to the geometrical problem such as the percolation even for d = 3. As mentioned in In- 
troduction, the transfer-matrix method has some advantages over the Monte-Carlo method. 
Actually, as for the d = 2 percolation, the transfer-matrix approach [17] has made a unique 
contribution, although its accuracy as to the critical indices is not particularly superior to 
that of Monte Carlo. Lastly, let us address a remark on the advantage of our approach: For 
example, according to Fortuin and Kasteleyn [lSj . the g-state- Potts model admits a "geo- 
metrical" representation in terms of the bond percolation. Namely, based on the percolation 
framework, one is able to extend the integral number q to a continuously variable one. In 
fact, in d = 3, an extensive Monte-Carlo study ^| reports an existence of the critical value 
q c = 2.45(10), above which the magnetic transition becomes discontinuous; see also Ref. j^J. 
However, the nature of this singularity is not fully understood, because the Monte Carlo 
method cannot deal with the complex-g number (due to the negative sign problem). On the 
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contrary, the transfer-matrix method is free this difficulty. Actually, with the transfer-matrix 
method, in d = 2, in the complex-g plane, the distribution of zeros (of the partition function) 



was investigated 



21 
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231 ] . As for d = 3, however, similar consideration has not yet been 



done due to lack of efficient algorithm. Our scheme meets such a requirement. Moreover, 



2J, and 



the low-lying spectrum of the Potts model in d = 3 is of fundamental significance 
the problem is also remained unsolved. In this sense, our scheme provides a step toward a 
series of such longstanding problems, to which the Monte-Carlo method does not apply. An 
effort toward this direction is in progress, and it will be addressed in future study. 
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(b) 



FIG. 1: (a) A drawing of the transfer-matrix bar. A unit of the transfer-matrix slice consists of N 
sites, (b) An example of connectivity among the N(= 4) sites. Such a connectivity is represented 
by the matrix notation, Eq. 
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FIG. 2: The scaled correlation length £/L is plotted for the percolation probability p with the 
fixed anisotropy parameter R = 3.3. The symbols +, x, *, □, ■. o, and • denote the system 
sizes N = 4, 5, 6, 7, 8, 9, and 10, respectively; note that the relation L = y/~N holds. From 
the intersection point of these curves, we read off the location of the critical point (percolation 
threshold) as p c w 0.096. 
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FIG. 3: The approximate critical point p c (L\,L2) is plotted for (2/(Li + L2)) 2 and the various 
values of the anisotropy parameter R = 2.8, 3.3, and 3.8. The symbols +, x, and * denote the 
differences of the system sizes N\ — N2 = 1, 2, and 3, respectively. We also presented the slopes 
with the least-square fit to the data as the dashed lines. We see that the finite-size corrections to 
p c are suppressed on setting R = 3.3. 
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FIG. 4: The approximate correlation-length critical exponent v(Li, L<i) is plotted for (2/(Li+L2)) 2 
with the fixed anisotropy parameter R = 3.3. The symbols +, x, and * denote the differences 
of the system sizes N± — N2 = 1, 2, and 3, respectively. The least-square fit to these data yields 
v = 0.812(15) in the limit L -> 00. 
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FIG. 5: The approximate beta function j3(L\,L2) (|19|) is plotted for the percolation probability 
p with the fixed anisotropy parameter R = 3.3. The symbols +, x, and * show that the pairs of 
the system sizes are (N±,N2) = (6,8), (7,9), and (8,10), respectively. We also presented a slope 
— (p — Pcjjv with p c = 0.0958 and v = 0.812 determined in Figs. Eland0J 
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FIG. 6: The approximate beta function f3(L\,L2) ifHfy is plotted for the percolation probability 
p at R = 1 (isotropic point). The symbols +, x, and * show that the pairs of the system sizes 
are (iVi, iVb) = (6,8), (7,9), and (8,10), respectively. We also presented a slope — {p — p c )/v with 
p c = 0.2488126 and v = 0.89 reported in Ref. dashed line. 
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